setwd("...")
options("modelsummary_format_numeric_latex" = "plain")
library(modelsummary)
library(kableExtra)
library(ggplot2)
library(tidyverse)
library(broom)
library(paletteer)
source("extra/func.R")
theme_set(my_theme())

data = read.csv("data/dataset.csv") %>%
  filter(ccaa == "catalunya")

### ====================================================================

## MAPS & CIA for TABLES

coef_recode = c(
  "patr33_bin" = "Employers' association",
  'CNT_bin' = 'CNT union',
  'UGT_bin' = 'UGT union',
  'pop1930_l' = 'Population (log)',
  'izq1936' = 'Leftist support 1936',
  'dias_control_rep_l' = 'Days of Republican control (log)',
  'analfab30_total' = 'Illiteracy 1930',
  'muerte_all_l' = 'Total Republican killings (log)',
  'clerics_all_l' = 'Number of clerics (log)',
  'patr33_bin:CNT_bin' = 'Employer $\\times$ CNT union',
  'patr33_bin:UGT_bin' = 'Employer $\\times$ UGT union'
)

gs = list(
  list("raw" = "nobs", "clean" = "$n$", "fmt" = 0),
  list("raw" = "r.squared", "clean" = "$R^2$", "fmt" = 2),
  list("raw" = "adj.r.squared", "clean" = "Adj. $R^2$", "fmt" = 2))

n = "+ p $<$ 0.1, * p $<$ 0.05, ** p $<$ 0.01, *** p $<$ 0.001. Outcome variable is the logged number of deaths for each professional category. Individuals working in law enforcement are excluded. Only municipalities in Catalonia."

### ====================================================================
## VIOLENCE BY OCCUPATION IN CATALUNYA

# New variables
data = data %>%
  mutate(
    muerte_leader_l = log(muerte_leader + 1),
    muerte_secular_l = log(muerte_clero_s + 1),
    muerte_regular_l = log(muerte_clero_r + 1),
    muerte_clero_l = log(muerte_clero_r + muerte_clero_s + 1),
    muerte_rest_l = log(muerte_other + 1),
    muerte_all_l = log(muerte_leader + muerte_clero_r + muerte_clero_s + muerte_other + 1),
    muerte_all_full_l = log(muerte_leader + muerte_clero_r + muerte_clero_s +
      muerte_security + muerte_other + 1))

data = data %>%
  mutate(muerte_all = muerte_leader + muerte_clero_r + muerte_clero_s + muerte_other) %>%
  mutate(share_leader = muerte_leader / muerte_all) %>%
  mutate(share_clero = (muerte_clero_r + muerte_clero_s) / muerte_all) %>%
  mutate(share_rest = muerte_other / muerte_all)


# Formula components
dv = c("muerte_clero_l", "muerte_leader_l", "muerte_rest_l")
ivars = c("patr33_bin", "CNT_bin", "UGT_bin")
controls = c("izq1936", "dias_control_rep_l", "analfab30_total",
  "pop1930_l", "clerics_all_l", "muerte_all_l")
form_right = paste0(paste(c(ivars, controls), collapse = " + "))

# Models
lm_cat = modnames(lapply(dv, function(x)
  lm(as.formula(paste0(x, " ~ ", form_right)), data = data)))

# Header for table
header = rep(1, length(lm_cat)+1)
names(header) = c(" ", gsub("muerte_(.*)_l", "\\1", dv) %>%
  recode("clero" = "Clergy", "rest" = "Rest", "leader" = "Local leaders"))

## Table
modelsummary(lm_cat,
  estimate = "{estimate}{stars}",
  output = "latex",
  coef_rename = coef_recode,
  coef_omit = "prov_name",
  gof_map = gs,
  title = "Linear models on Republican victimization in Catalonia by professional category\\label{tab:lm_cat_profesiones}",
  threeparttable = TRUE,
  escape = FALSE) %>%
add_header_above(header) %>%
kable_styling(latex_options = c("hold_position")) %>%
footnote(general = n, threeparttable = TRUE, footnote_as_chunk = TRUE, escape = FALSE) %>%
save_kable(file = "output/tab_lm_catalunya_categorias.tex")

## Plot
# Get coefficient estimates
coef_df = modelplot(lm_cat, draw = FALSE, coef_map = c(
    "patr33_bin" = "Employers' association",
    "CNT_bin" = "CNT union", "UGT_bin" = "UGT union")) %>%
  mutate(model_lab = recode(model, "(1)" = dv[1], "(2)" = dv[2], "(3)" = dv[3])) %>%
  mutate(model_lab = recode(model_lab, "muerte_clero_l" = "Clergy",
    "muerte_leader_l" = "Local leaders", "muerte_rest_l" = "Rest of civilians")) %>%
  rename(est = estimate, se = std.error) %>%
  mutate(
    upr95 = est + qnorm(0.975) * se,
    upr90 = est + qnorm(0.95) * se,
    lwr95 = est - qnorm(0.975) * se,
    lwr90 = est - qnorm(0.95) * se)

p = ggplot(coef_df %>% filter(grepl("Employers", term)),
    aes(x = model_lab, y = est)) +
  geom_point(position = position_dodge(1/5)) +
    geom_hline(yintercept = 0, linetype = "dotted") +
  geom_errorbar(aes(ymin = lwr90, ymax = upr90), width = 0, linewidth = 1.1,
    position = position_dodge(1/5)) +
  geom_errorbar(aes(ymin = lwr95, ymax = upr95), width = 0, position = position_dodge(1/5)) +
  labs(x = "\nOutcome", y = "Est. of Employers' associations\n(and 90/95% CIs)\n") +
  theme(legend.title = element_blank(), legend.position = "bottom") +
  scale_color_paletteer_d(`"wesanderson::Darjeeling1"`)
ggsave("output/coefplot_categorias_catalunya.pdf",
  height = 3, width = 4.5, device = "pdf")
